Pascal Super Library/Pascal Super Library (CW International)(1997).bin/MATH/PASCMP/PASFFT.PAS

Pascal/Delphi Source File | 1994-06-29 | 5KB | 211 lines

{ Complex Fast Fourier Transform Unit for Borland Pascal      }
{ (as an example for using complex mathematics unit PasCmplx) }
{ (c)1994 by Alex Klimovitski                                 }

{ Based on algorithm from:                                    }
{ Marple, S.L. Digital Spectral Analysis. Prentice-Hall, Inc.,}
{   Englewood Cliffs, N.J., 1988.                             }

{ See FFTTEST.PAS for examples of using.                      }

unit PasFFT;

{$N+}
{$G+}

interface

uses PasCmplx;

type
  PCmplxArr = ^TCmplxArr;
  TCmplxArr = array[1..4096] of Complex;

function CFFTInit(N: Integer; IsInverse: Boolean): Boolean;
{initializes the unit for the transform of N points.
 N -- number of points. N MUST be an integer power of 2,
   for ex.: 64, 1024, 4096.
 IsInverse -- false for direct Fourier tramsform,
   true for inversed Fourier transform.
 The function returns true if initialization was successful.
 NOTE: You don't need to call this function explicitly.
}

procedure CFFTClear;
{returns the memory occupied by internal tables to the heap.
 NOTE: You can call this function to free memory if you don't
   need to perform more transforms. This function is called
   automatically when the program ends.
}

function CFFT(N: Integer; IsInverse: Boolean; X: PCmplxArr; T: Double): Boolean;
{performs Fast Fourier Transform.
 N -- number of points. N MUST be an integer power of 2.
 IsInverse -- false for direct Fourier tramsform,
   true for inversed Fourier transform.
 X -- Pointer on an array of N complex numbers to be transformed.
 T -- time step
 The result is placed in array pointed by X.
 The function returns true if transform was successful and false
   if it failed.
 NOTE: Subsequent transforms with the same number of points N and
   in the same direction IsInverse are performed quicker because
   of using of the same internal table.
}

implementation

const
  GW: PCmplxArr = nil;

var
  GN: Integer;
  GNExp: Integer;
  GIsInverse: Boolean;

function CFFTInit(N: Integer; IsInverse: Boolean): Boolean;
var
  P, Q: Complex;
  S, C: Double;
  I: Integer;
begin
  CFFTInit := false;

  GNExp := 0;
  I := 1;
  while I < N do
  begin
    I := I * 2;
    Inc(GNExp);
  end;
  if I <> N then Exit;

  GN := N;
  GIsInverse := IsInverse;

  S := 8.0 * ArcTan(1.0) / GN;
  CSinCosR(S, S, C);
  P := Cmplx(C, S);

  if IsInverse then P := Conjug(P);

  GetMem(GW, GN * SizeOf(Complex));
  if GW = nil then Exit;
  Q := C1;

  for I := 1 to GN do
  begin
    GW^[I] := Q;
    Q := CMul(Q, P);
  end;

  CFFTInit := true;
end;

procedure CFFTClear;
begin
  if GW <> nil then
    FreeMem(GW, GN * SizeOf(Complex));
  GN := 0;
end;

function CFFT(N: Integer; IsInverse: Boolean; X: PCmplxArr; T: Double): Boolean;
var
  P, Q: Complex;
  I, J, K, JJ, KK, LL, MM, NN: Integer;
begin
  CFFT := false;

  if X = nil then Exit;

  if (GW = nil) or (GN <> N) or (GIsInverse <> IsInverse) then
  begin
    CFFTClear;
    if not CFFTInit(N, IsInverse) then Exit;
  end;

  MM := 1;
  LL := GN;

  for K := 1 to GNExp do
  begin
    NN := LL div 2;
    JJ := MM + 1;

    I := 1;
    while I <= GN do
    begin
      KK := I + NN;
      P := CAdd(X^[I], X^[KK]);
      X^[KK] := CSub(X^[I], X^[KK]);
      X^[I] := P;
      Inc(I, LL);
    end; {while I}

    {if NN = 1 then Continue;}

    for J := 2 to NN do
    begin
      Q := GW^[JJ];

      I := J;
      while I <= GN do
      begin
        KK := I + NN;
        P := CAdd(X^[I], X^[KK]);
        X^[KK] := CMul(CSub(X^[I], X^[KK]), Q);
        X^[I] := P;
        Inc(I, LL);
      end; {while I}
      JJ := JJ + MM;
    end; {for J}

    LL := NN;
    MM := MM * 2;
  end; {for K}

  NN := GN div 2;
  MM := N - 1;
  J := 1;

  for I := 1 to MM do
  begin
    if I < J then
    begin
      P := X^[J];
      X^[J] := X^[I];
      X^[I] := P;
    end; {if I < J}

    K := NN;
    while K < J do
    begin
      Dec(J, K);
      K := K div 2;
    end; {while K < J}

    Inc(J, K);
  end; {for I}

  if IsInverse then P := Cmplx(1.0 / (T * GN), 0.0)
  else P := Cmplx(T, 0.0);

  for I := 1 to GN do
    X^[I] := CMul(X^[I], P);

  CFFT := true;
end;

var
  ExitSave: Pointer;

procedure CFFTExit; far;
begin
  ExitProc := ExitSave;
  CFFTClear;
end;

begin
  ExitSave := ExitProc;
  ExitProc := @CFFTExit;
end.